%CFK_eDV   is a modified version of the original CFK model [1] for curve-fitting and estimation of DLCO and VCO, using a height-weight-sex-based estimate of blood volume

% Copyright (c) 2020, Tilo Winkler and R. Scott Harris
%
%   This program is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
%
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
%
%   You should have received a copy of the GNU General Public License
%   along with this program.  If not, see <https://www.gnu.org/licenses/>.
%
%
%
% ____REFERENCES____
% 1. Peterson JE, Stewart RD. Predicting the carboxyhemoglobin levels resulting from carbon monoxide exposures. Journal of Applied Physiology. 1975 Oct 1;39(4):633–8
% 2. Nadler SB, Hidalgo JH, Bloch T. Prediction of blood volume in normal human adults. Surgery. 1962 Feb;51(2):224-32.
%      Nadler Equation[3]:
%          Men: Blood Volume   = (0.3669 × H^3) + (0.03219 × W) + 0.6041
%          Women: Blood Volume = (0.3561 × H^3) + (0.03308 × W) + 0.1833
%      source: Physiology, Blood Volume - StatPearls - NCBI Bookshelf -- https://www.ncbi.nlm.nih.gov/books/NBK526077/
%
%
% ____INPUT PARAMETERS____
%
%  Name         Meaning (unit)                      Default
%
%  in.ID        Subject ID
%  in.Weight    (kg)                                 75
%  in.Height    (m)
%  in.Sex       (1:male; 0:female)
%  in.FIO2      Inspiratory oxygen fraction (0..1)    0.21
%  in.MV        Minute ventilation (l/min):           5
%  in.VDVT      VD/VT ratio                           0.25
%  in.ctHb      Hb (g/dL)                            11.6
%  in.HbCO_0    Baseline HbCO (%)                     0         -- obsolete for DV model
%  in.COppm     Inhaled CO Conc. (ppm)              200
%  in.PB        Barometric pressure (mmHg)          760
%  in.HbCO_m    Measured HbCO fraction (%) at time(s) t_m (see below)
%  in.t_m       Time of measurement relative to start of CO (min)


function res = CFK_eDV( uin, ts, model_type)

% ____INIT____
male = 1;
VdotCO_0 = 0.007;   % CO production set at 0.007 mL/min
DLCO_0 = 5;   % start value for estimation


% ____MAIN____
res = nCalc;


% ____NESTED FUNCTIONS____
% NOTE NAMING CONVENTION: 'n'FunctionName for nested functions, e.g., nCalc or nHbCO_residuals sharing variables with the parent function


function res = nCalc % -----------------------------------------------
%nCalc   estimate model parameters

global HbCO_0_est HbCO_est residuals % for additional return values from  .. = fminsearch( @nHbCO_residuals, ..)

if uin.Sex == male, % height, weight, and sex-based blood volume estimate [2]
    uin.Vb = 1000 * ( (0.3669 * uin.Height^3) + (0.03219 * uin.Weight) + 0.6041 ); % Vb in ml
  else
    uin.Vb = 1000 * ( (0.3561 * uin.Height^3) + (0.03308 * uin.Weight) + 0.1833 ); % Vb in ml
    end

uin.VdotA = uin.MV * ( 1 - uin.VDVT) * 1000;   % alveolar ventilation in ml/min
uin.PICO = ( uin.COppm * ( uin.PB - 47)) / 1000000;   % convert inhaled CO from ppm to partial pressure (mmHg)
uin.cHbO2max = 1.38 * uin.ctHb / 100;   % maximum O2 (or CO) in blood (mL O2 or CO/mL blood)
uin.cHbCO_0 = uin.HbCO_0 / 100 * uin.cHbO2max;   % cHbCO at the start, in mL CO/mL blood
uin.cHbCO_m = uin.HbCO_m / 100 * uin.cHbO2max;   % cHbCO at time t_m, in mL CO/mL blood

switch model_type
  case 'DV'
    p_0 = [ DLCO_0 VdotCO_0];
  case 'D'
    p_0 = [ DLCO_0];
  case 'V'
    p_0 = [ VdotCO_0];
  otherwise
    error('Unknown model type');
    end

HbCO_est = NaN( size( uin.t_m_all));
residuals = HbCO_est;
opt = optimset( 'MaxFunEvals', 100000, 'TolX', 1e-7, 'TolFun', 1e-7);
p_est = fminsearch( @nHbCO_residuals, p_0, opt); % estimate DLCO

switch model_type
  case 'DV'
    res.DLCO_est   = p_est(1); % second position in X vector !
    res.VdotCO_est = p_est(2);
  case 'D'
    res.DLCO_est   = p_est(1); % second position in X vector !
    res.VdotCO_est = VdotCO_0;
  case 'V'
    res.DLCO_est   = uin.predicted_DLCO; % second position in X vector !
    res.VdotCO_est = p_est(1);
    end
res.HbCO_0_est = HbCO_0_est;
res.HbCO_est = HbCO_est;
res.residuals = residuals;
    
end % -- of function --------------------------------------


function err = nHbCO_residuals( p_est) % ---------------------------------
%nHbCO_residuals( p_est)   calculate residuals (predicted - measured HbCO) for the estimated parameters p_est

global HbCO_0_est HbCO_est residuals

switch model_type
  case 'DV'
    DLCO_est    = p_est(1);
    VdotCO_est  = p_est(2);
  case 'D'
    DLCO_est    = p_est(1);
    VdotCO_est  = VdotCO_0;
  case 'V'
    DLCO_est    = uin.predicted_DLCO;
    VdotCO_est  = p_est(1);
  otherwise
    error('Unknown model type');
    end

% constrains for failed/non-converging estimates
if ( DLCO_est < 0.0001) || ( DLCO_est > 1000), % relax constrains 
    err = 10;
    return
    end
if ( VdotCO_est < 0.0001) || ( VdotCO_est > 1), % relax constrains
    err = 10;
    return
    end
    
[ cHbCO_0_est, HbCO_0_est] = nCFK_cHbCO_0_est( DLCO_est, VdotCO_est);

for k = 1:length( uin.t_m_all)
    HbCO_est(k) = nCFK_FiO2( cHbCO_0_est, DLCO_est, VdotCO_est, uin.t_m_all(k)); end % estimated content
residuals = HbCO_est - uin.HbCO_m_all;
err = sum( residuals.^2); % least square error of estimates compared to measurements
end % n-function: nHbCO_residuals


function [ cHbCO_0_est, HbCO_0_est] = nCFK_cHbCO_0_est( DLCO_i, VdotCO_i) % -------------------------------
%nCFK_cHbCO_0_est    estimate baseline carboxyhemoglobin content cHbCO and HbCO assuming a steady state for the given DLCO and VdotCO prior to iCO (--> PICO = 0)

cHbCOtold = 0.0;
cHbO2 = uin.cHbO2max - cHbCOtold;

d = 1.0; % absolute error
while d >= 0.00001 % difference between old and new values of HbCO, stop when difference is small
    c = B( DLCO_i, uin.PB, uin.VdotA) * VdotCO_i;
    cHbCOtnew = c / A( cHbO2, uin.PB, uin.FIO2, 0);
    d = abs( cHbCOtnew - cHbCOtold);
    cHbCOtold = cHbCOtnew;
    cHbO2 = uin.cHbO2max * ( 1 - cHbCOtold / ( cHbCOtold + cHbO2));
    end

cHbCO_0_est =  cHbCOtnew;
HbCO_0_est = ( cHbCOtnew * 100) / uin.cHbO2max;
end % n-function: nCFK_cHbCO_0_est


function HbCO_est = nCFK_FiO2( cHbCO_0_i, DLCO_i, VdotCO_i, t_i) % -------------------------------
%nCFK_FIO2    calculate estimate of carboxyhemoglobin

cHbCOtold = 0.0;
cHbO2 = uin.cHbO2max - cHbCOtold;

d = 1.0; % absolute error
while d >= 0.00001 % difference between old and new values of HbCO, stop when difference is small
    p = ( -t_i * A( cHbO2, uin.PB, uin.FIO2, uin.PICO)) / ( uin.Vb * B( DLCO_i, uin.PB, uin.VdotA));
    a = exp( p);
    b = A( cHbO2, uin.PB, uin.FIO2, uin.PICO) * cHbCO_0_i   - B( DLCO_i, uin.PB, uin.VdotA) * VdotCO_i - uin.PICO;
    c = B( DLCO_i, uin.PB, uin.VdotA) * VdotCO_i + uin.PICO;
    cHbCOtnew = ( a * b + c) / A( cHbO2, uin.PB, uin.FIO2, uin.PICO);
    d = abs( cHbCOtnew - cHbCOtold);
    cHbCOtold = cHbCOtnew;
    cHbO2 = uin.cHbO2max * ( 1 - cHbCOtold / ( cHbCOtold + cHbO2));
    end

HbCO_est = ( cHbCOtnew * 100) / uin.cHbO2max;
end % n-function: nCFK_FiO2

end % MAIN FUNCTION ------------------------------------------------------------


% ____SUB-FUNCTIONS____

function [A] = A( cHbO2, PB, FIO2, PICO) % -------------------------------------
%A   calculate subfunction A
 
M = 218; % ratio of the affinity of blood for CO to that for O2
PIO2 = FIO2 * ( PB - 47 - PICO);
if PIO2 >= 149; 
    PcO2 = PIO2 - 49;               
  else
    PcO2 = 1 / ( 0.072 - 0.00079 * PIO2 + 0.000002515 * PIO2 ^2);
    end
A = PcO2 / ( M * cHbO2);
end % s-function


function [B] = B( DLCO, PB, VdotA) % ------------------------------------------
%B   calculate subfunction B
 
PL = PB - 47;
B = ( 1 / DLCO) + ( PL / VdotA);
end % s-function
